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ABSTRACT 


The  specific  objective  of  this  research  is  the  development  and  improvement  of 
theoretical  models  to  simulate  the  effect  of  nonsphericity  on  single-scattering  properties  of 
cirrus  cloud  particles  in  the  visible  and  infrared  spectral  regions.  First,  we  have  shown  that 
using  a  matrix  inversion  scheme  based  on  a  special  LU  factorization  rather  than  on  the 
standard  Gaussian  elimination  significantly  improves  the  numerical  stability  of  T-matrix 
computations  for  nonabsorbing  and  weakly  absorbing  nonspherical  particles.  As  a  result, 
the  maximum  convergent  size  parameter  for  particles  with  small  or  zero  absorption  can 
increase  by  a  factor  of  several  and  can  exceed  100.  Comparisons  of  T-matrix  and  geometric 
optics  (GO)  computations  for  large,  randomly  oriented  spheroids  and  finite  circular  cylinders 
show  that  the  range  of  applicability  of  the  ray  tracing  approximation  depends  on  the  • 
imaginary  part  of  the  refractive  index  and  is  different  for  different  elements  of  the  scattering 
matrix. 

Second,  we  use  exact  T-matrix  computations  and  physical  considerations  based  on 
the  Kirchhoff  approximation  to  show  that  the  <f-function  transmission  peak  predicted  by  the 
GO  approximation  for  hexagonal  ice  crystals  is  an  artifact  of  GO  completely  ignoring 
physical  optics  effects  and  must  be  convolved  with  the  Fraunhofer  pattern,  thereby 
producing  a  phase  function  component  with  an  angular  profile  similar  to  the  standard 
diffraction  component.  We  have  performed  this  convolution  with  a  simple  procedure  which 
supplements  the  standard  ray  tracing  code  and  makes  the  computation  of  the  phase  function 
and  its  Legendre  expansion  both  more  physically  realistic  and  more  accurate. 
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INTRODUCTION 


This  research  project  is  intended  to  deliver  crucial  modeling  products  for  use  in  the 
project  headed  by  Dr.  B.-C.  Gao  and  called  "Cirrus  Cloud  Characterization  and  Correction." 

It  is  well  known  that  the  accuracy  of  determining  the  cirrus  cloud  optical  thickness  and 
correcting  the  effect  of  cirrus  on  the  retrieved  aerosol  and  lower-cloud  optical  thicknesses 
and  surface  reflectance  critically  depends  on  the  accuracy  of  modeling  the  cirrus  cloud 
particle  single-scattering  properties,  first  of  all  the  scattering  phase  function.  The  knowledge 
of  the  phase  function  is  also  crucially  important  in  an  accurate  evaluation  of  the  radiative 
effect  of  cirrus  on  climate.  Cirrus  clouds  consist  of  ice  particles  having  an  wide  variety  of 
shapes.  For  typical  cirrus  environments  ice  crystals  exist  both  in  single  and  aggregated 
forms,  and  the  surface  of  the  ice  particles  can  be  highly  irregular.  Because  of  ignorance  • 
about  the  distribution  of  cirrus  particles  over  shapes,  knowledge  of  cirrus  scattering  phase 
function  has  progressed  very  slowly.  Theoretical  phase  function  calculations  have  mostly 
assumed  that  cirrus  ice  particles  are  spheres  or  randomly-oriented,  infinitely-long  circular 
cylinders;  this  remains  a  common  assumption  both  in  global  climate  models  and  in  satellite 
retrieval  algorithms.  The  geometric-optics  hexagonal-crystal  phase  functions  of  Takano  and 
Liou  have  become  quite  popular,  but  do  not  seem  to  fit  the  known  facts  very  well.  Phase 
function  calculations  for  small  ice  particles  have  been  difficult  and  computationally 
intensive. 

Because  of  these  factors,  the  focus  of  this  research  during  the  period  23  September 
1996  -  30  September  1997  has  been  on  further  development  and  use  state-of-the-art 
numerical  techniques  for  computing  scattering  properties  of  nonspherical  ice  crystals.  The 
specific  objective  has  been  the  development  and  improvement  of  theoretical  models  to 
simulate  the  effect  of  nonsphericity  on  single-scattering  properties  of  cirrus  cloud  particles 
in  the  visible  and  infrared  spectral  regions. 
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IMPROVEMENTS  OF  THE  T-MATRIX  METHOD 


The  Waterman's  T-matrix  approach  is  one  of  the  most  powerful  exact  techniques  for 
computing  light  scattering  by  nonspherical  particles  based  on  solving  Maxwell's  equations. 
However,  even  this  method  can  exhibit  convergence  problems  when  any  of  the  variables 
defining  the  scattering  particle  (size  parameter,  deviation  from  sphericity,  refractive  index) 
becomes  too  extreme.  Standard  T-matrix  computations  become  especially  ill-conditioned 
for  particles  with  a  small  or^zero  imaginary  part  of  the  refractive  index  because  of  the  strong 
effect  of  the  ripple  structure.  For  example,  the  maximum  convergent  equivalent-sphere  size 
parameter  (i.e.,  the  ratio  of  particle  circumference  to  wavelength  of  scattered  light  for  the 
surface-equivalent  sphere)  in  double-precision  FORTRAN  computations  for  oblate  spheroids 
with  refractive  index  1.53  and  aspect  ratio  2  is  only  8,  whereas  the  maximum  convergent  * 
size  parameter  for  the  same  particles  but  with  a  moderately  absorbing  refractive  index  of 
1.53  +  0.001  i  is  33. 

The  origin  of  the  numerical  instability  of  the  standard  T-matrix  procedure  for  extreme 
values  of  particle  characteristics  can  be  explained  as  follows.  Calculations  based  on  the 
extended  boundary  condition  method  (EBCM)  assume  the  representation  of  the  T-matrix  in 
the  form  7  =  -Q~^  RgQ,  where  the  elements  of  the  matrices  Q  and  RgQ  are  integrals  over 
the  particle  surface.  The  numerical  inversion  of  the  matrix  Q  is  usually  performed  using  the 
standard  Gaussian  elimination  (GE).  Unfortunately,  the  calculation  of  the  inverse  matrix 
Q"’  is  an  ill-conditioned  procedure  strongly  affected  by  round-off  errors  and  by  the  fact  that 
different  elements  of  the  Q  matrix  can  differ  by  many  orders  of  magnitude.  As  a  result,  T- 
matrix  computations  for  extreme  particle  parameters  can  be  poorly  convergent  and  even 
divergent. 

This  sensitivity  of  the  standard  T-matrix  procedure  to  weak  or  zero  absorption  has 
been  a  serious  limiting  factor  since  many  commonly  encountered  substances  are  weakly 
absorbing  in  some  spectral  ranges.  The  example  of  water  ice  in  the  visible  is  particularly 
important  because  quite  often  a  significant  fraction  of  cirrus  and  contrail  ice  particles  are  not 
much  larger  than  a  visible  wavelength,  thus  potentially  precluding  the  use  of  the  geometric 
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optics  approximation  (GO)  in  light  scattering  computations. 

During  the  course  of  this  research,  we  have  significantly  improved  the  T-matrix 
technique  for  computing  light  scattering  by  nonspherical  ice  crystals.  Specifically,  we  have 
shown  that  using  a  matrix  inversion  scheme  based  on  a  special  LU  factorization  rather  than 
on  the  standard  Gaussian  elimination  significantly  improves  the  numerical  stability  of  T- 
matrix  computations  for  nonabsorbing  and  weakly  absorbing  nonspherical  particles.  We 
have  also  developed  an  improved  scheme  for  evaluating  Clebsch-Gordon  coefficients  with 
large  quantum  numbers  which  allowed  us  to  extend  the  analytical  orientational  averaging 
method  developed  by  Mishchenko  D-  Opt.  Soc.  Am.  A  8,  871  (1991)]  to  larger  size 
parameters.  As  a  result,  the  maximum  convergent  size  parameter  for  particles  with  small 
or  zero  absorption  can  increase  by  a  factor  of  several  and  can  exceed  100.  The  ability  of 
the  new  T-matrix  scheme  to  treat  such  large  size  parameters  is  accompanied  by  an  extremely  . 
high  numerical  efficiency  which  makes  our  code  orders  of  magnitude  faster  than  any 
alternative  technique  for  exactly  computing  light  scattering  by  nonspherical  particles  in 
random  orientation.  The  unique  capabilities  of  the  T-matrix  code  make  it  very  useful  in 
practice,  but  also  make  difficult  checks  of  its  numerical  accuracy  since  independent  results 
for  the  largest  size  parameters  cannot  be  obtained  with  any  other  currently  available 
method.  Therefore,  we  have  paid  special  attention  to  making  sure  that  our  calculations  fully 
satisfy  such  fundamental  physical  constraints  as  symmetry,  reciprocity,  and  energy 
conservation.  Also,  we  have  computed  benchmark  results  for  a  challenging  test  case  that 
can  be  used  for  checking  the  accuracy  of  the  most  advanced  nonspherical  scattering  codes 
at  higher  frequencies. 

Using  the  present  version  of  the  T-matrix  code,  we  have  extended  comparisons  of 
exact  T-matrix  and  approximate  ray  tracing  calculations  to  much  larger  size  parameters  and 
to  all  elements  of  the  scattering  matrix.  Our  results  suggest  that  equivalent-sphere  size 
parameters  larger  than  about  80  are  already  big  enough  to  ensure  acceptable  accuracy  of 
GO  phase  function  computations  (except,  perhaps,  at  exactly  the  backscattering  direction). 
However,  GO  calculations  of  the  other  elements  of  the  scattering  matrix  are  stronger 
affected  by  wave  effects  and  become  reasonably  accurate  only  at  significantly  larger  size 
parameters.  GO  calculations  of  lidar  depolarization  can  be  expected  to  be  especially 
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inaccurate  unless  the  equivalent-sphere  size  parameter  exceeds  several  hundred. 

An  interesting  result  of  our  calculations  is  that  the  T-matrix  method  can  be 
successfully  applied  to  large  sharp-edged  particles  such  as  finite  circular  cylinders  with 
equivalent-sphere  size  parameters  exceeding  100.  It  has  been  often  claimed  that  the 
presence  of  sharp  edges  can  be  difficult  to  handle  with  a  method  that  uses  "smooth" 
spherical  functions  in  the  internal  and  scattered  field  decompositions.  We  have  found, 
however,  that  the  use  of  a  special  numerical  integration  scheme  for  computing  the  surface 
integrals  needed  to  calculate  the  T  matrix  ameliorates  the  problem  of  sharp  edges  and  makes 
T-matrix  computations  for  cylinders  almost  as  accurate  as  those  for  surface-  and  aspect-ratio- 
equivalent  smooth-shaped  spheroids.  Much  more  difficult  problems  are  encountered  when 
the  T-matrix  method  is  applied  to  particles  with  large  aspect  ratios.  In  this  case  a  single 
spherical  function  expansion  of  the  internal  and  scattered  fields  can  fail,  and  the  use  of  . 
several  overlapping  subdomain  spherical  function  expansions  may  become  necessary. 

This  research  was  published  in  Applied  Optics  (see  references  below). 


INCORPORATION  OF  PHYSICAL  OPTICS  EFFECTS  INTO 
RAY-TRACING  PHASE  FUNCTIONS  COMPUTED  FOR  HEXAGONAL  ICE 

CRYSTALS 


It  is  well  known  that  a  convenient  way  of  representing  the  scattering  phase  function 
P(0)  for  aerosol  and  cloud  particles  is  expanding  it  in  Legendre  polynomials  as 


P(0) 


max 

E  Pn(COS0), 

n=0 


(1) 


where  ©  is  the  scattering  angle,  P„(cos  0)  are  Legendre  polynomials,  and  the  value  of  the 
upper  summation  limit  depends  on  the  desired  numerical  accuracy  of  the  expansion. 
Since  the  number  of  numerically  significant  terms  in  the  Legendre  expansion  is  finite  and 
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often  relatively  small,  this  expansion  can  be  used  for  efficiently  computing  the  phase 
function  .for  essentially  any  number  of  scattering  angles  with  a  small  consumption  of  CPU 
time.  Furthermore,  the  Legendre  expansion  coefficients  can  be  used  to  directly  compute 
the  Fourier  components  of  the  phase  function  via  simple  and  exact  analytical  formulas, 
which  is  the  first  step  in  radiative  transfer  computations  using  different  numerical  techniques. 

The  Legendre  expansion  coefficients  for  the  widely  used  Flenyey-Greenstein 
phase  function  are  given  by  the  simple  analytical  formula 

=  {2n  +  ^)g^,  (2) 

where  g  is  the  asymmetry  parameter.  Efficient  analytical  methods  based  on  solving 
Maxwell's  equations  exist  for  computing  the  expansion  coefficients  for  spherical  particles 
and  randomly  oriented,  rotationally  symmetric  nonspherical  particles.  For  irregular  particles 
with  sizes  much  larger  than  the  wavelength  of  the  incident  radiation,  such  as  cirrus  cloud 
particles  in  the  visible,  direct  numerical  solutions  of  Maxwell's  equations  do  not  currently 
exist.  Therefore  the  expansion  coefficients  have  to  be  computed  using  an  approximate 
technique  such  as  the  geometric  optics  approximation  (GO).  Using  the  orthogonality 
property  of  Legendre  polynomials,  we  easily  derive  from  equation  (1) 

TT 

x„  =  f  d0P(0)P„(co50)sin0.  (3) 

The  integral  in  equation  (3)  can  be  calculated  numerically  by  using  a  quadrature  formula 
provided  that  the  phase  function  values  at  the  division  points  are  known.  This  numerical 
approach  works  well  if  the  phase  function  is  rather  smooth  but  becomes  problematic  for 
particles  having  parallel  planes  such  as  hexagonal  columns  and  plates  or  finite  circular 
cylinders.  In  this  case,  the  standard  GO  predicts  a  strong,  infinitesimally  narrow  peak  in  the 
exact  forward-scattering  direction  which  is  caused  by  rays  that  undergo  two  refractions 
through  parallel  plane  facets  and  is  superimposed  on  the  diffraction  component  of  the  phase 
function.  This  effect  was  called  by  Takano  and  Liou  the  ^-function  transmission. 
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It  is  obvious,  however,  that  GO  predicts  the  infinitesimally  narrow  (^-function 
transmission  peak  only  because  it  completely  ignores  physical  optics  effects.  Simple 
physical  optics  considerations  cause  us  to  conclude  that  although  a  strong  non-diffraction 
forward-scattering  peak  does  exist  and  can  be  qualitatively  explained  in  GO  terms  as  a 
manifestation  of  the  (5-function  transmission,  it  nonetheless  has  an  appreciable  angular  width 
comparable  to  that  of  the  Fraunhofer  diffraction  peak  and  a  diffraction-like  angular  profile. 

It  is  clear  that  however  large  a  particle  is  compared  to  the  wavelength,  physical  optics  effec^ts  ■ 
will  preclude  the  appearance  of  perfect  singularities  in  the  scattering  pattern  like  the  <5- 
function  transmission  peak  in  the  phase  function.  Instead,  a  wave  front  emerging  from  a  flat 
crystal  facet  should  spread  and  produce  an  angular  intensity  distribution  in  the  far-field  zone 
similar  to  the  well  known  Fraunhofer  diffraction  pattern.  Of  course,  in  the  theoretical  limit 
of  an  infinite  size  parameter  the  <5-function  transmission  peak  becomes  a  true  ^-function.  . 
However,  the  angular  width  of  the  ^-function  transmission  peak  is  always  comparable  to  that 
of  the  Kirchhoff  diffraction  component,  and  as  long  as  the  latter  is  computed  explicitly,  so 
should  be  the  <5-function  transmission  component. 

We  used  exact  T-matrix  computations  for  rather  large  nonspherical  particles  to 
demonstrate  that  the  effect  that  can  be  interpreted  in  geometric  optics  terms  as  ^-function 
transmission  through  parallel  planes  indeed  results  in  a  quasi-Fraunhofer  forward-scattering 
peak  rather  than  in  a  true  <5-function  peak.  We  then  developed  a  very  simple  numerical 
procedure  which  incorporates  this  physical  optics  effect  in  the  standard  ray  tracing 
procedure  for  computing  the  phase  function  for  large  particles  having  parallel  plane  facets. 
This  procedure  not  only  makes  ray  tracing  computations  more  physically  relevant,  but  also 
simplifies  and  makes  more  accurate  the  computation  of  the  phase  function  and  its  Legendre 
expansion.  Although  the  accuracy  of  our  procedure  cannot  be  assessed  directly  due  to  the 
lack  of  exact  theoretical  methods  based  on  solving  Maxwell's  equations  and  applicable  to 
size  parameters  exceeding  several  hundred,  our  approach  is  physically-based  and  appears 
to  be  simple  and  well  justified  since  it  consists  of  directly  computing  the  amount  of  energy 
contained  in  the  <5-function  transmission  peak  and  convolving  it  with  the  Fraunhofer  angular 
pattern. 

As  an  example,  in  Figure  1  we  show  the  total  phase  function  for  polydisperse. 
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randomly  oriented  hexagonal  columns  with  length-to-diameter  ratio  2  and  distribution  of 
surface-equivalent-sphere  radii  given  by  the  standard  power  law: 


n(r)  = 


2r  2 

^2  ^-3 


2  2 

^2-^1 


r,  <  r  < 


0  otherwise, 


(4) 


The  parameters  and  r2  are  chosen  such  that  the  effective  radius  and  effective  variance  of 
the  distribution,  are  r^ff  =  40  //m  and  =  0.2.  The  refractive  index  is  1.3082  + 
/0.1 328  X 1 0  and  the  wavelength  is  >1  =  0.645  //m.  The  solid  curve  shows  the  original 
phase  function,  while  the  dotted  curve  shows  the  result  of  evaluating  the  Legendre 
expansion  of  equation  (1)  with  =  1000  terms.  It  is  seen  that  the  original  phase  ‘ 
function  and  the  Legendre  expansion  almost  perfectly  coincide  (relative  differences  less  than 
10"^).  The  forward-scattering  value  for  this  phase  function  is  1.1 60x10^  and  the  asymmetry 
parameter  is  0.8209. 

This  research  has  been  summarized  in  a  paper  submitted  to  Journal  of  Geophysical 
Research  (Atmospheres). 
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Figure  1.  Scattering  phase  function  for  polydisperse,  randomly  oriented  hexagonal  ice 
columns.  The  solid  curve  shows  the  original  phase  function  and  the  dotted  curve  shows  the 
result  of  evaluating  the  Legendre  expansion  of  equation  (1)  with  1000  terms. 
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infrared  spectral  regions.  First,  we  have  shown  that  using  a  matrix  inversion  scheme  based  on  a  special  LU 
factorization  rather  than  on  the  standard  Gaussian  elimination  significantly  improves  the  numerical  stability  of 
T-matrix  computations  for  nonabsorbing  and  weakly  absorbing  nonspherical  particles.  As  a  result  the 
maximum  convergent  size  parameter  for  particles  with  small  or  zero  absorption  can  increase  by  a  factor  of 
several  and  can  exceed  1 00.  Comparisons  of  T-matrix  and  geometric  optics  (GO)  computations  for  large 
randomly  oriented  spheroids  and  finite  circular  cylinders  show  that  the  range  of  applicability  of  the  ray  tracing 
approximation  depends  on  the  imaginary  part  of  the  refractive  index  and  is  different  for  different  elements  of 
the  scattering  matrix.  Second,  we  use  exact  T-matrix  computations  and  physical  considerations  based  on  the 
Kirchhoff  approximation  to  show  that  the  <f-function  transmission  peak  predicted  by  the  GO  approximation  for 
I  hexagonal  ice  crystals  is  an  artifact  of  GO  completely  ignoring  physical  optics  effects  and  must  be  convolved 
;  with  the  Fraunhofer  pattern,  thereby  producing  a  phase  function  component  with  an  angular  profile  similar  to 
the  standard  diffraction  component.  We  have  performed  this  convolution  with  a  simple  procedure  which 
supplements  the  standard  ray  tracing  code  and  makes  the  computation  of  the  phase  function  and  its  Legendre 
expansion  both  more  physically  realistic  and  more  accurate. 
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